Connectivity percolation in suspensions of hard platelets 
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Abstract 

We present a study on connectivity percolation in suspensions of hard platelets by means of 
Monte Carlo simulation. We interpret our results using a contact-volume argument based on 
an effective single-particle cell model. It is commonly assumed that the percolation threshold 
of anisotropic objects scales as their inverse aspect ratio. While this rule has been shown to 
hold for rod-like particles, we find that for hard plate-like particles the percolation threshold is 
non-monotonic in the aspect ratio. It exhibits a shallow minimum at intermediate aspect ratios 
and then saturates to a constant value. This effect is caused by the isotropic-nematic transition 
pre-empting the percolation transition. Hence the common strategy to use highly anisotropic, 
conductive particles as fillers in composite materials in order to produce conduction at low filler 
concentration is expected to fail for plate-like fillers such as graphene and graphite nanoplatelets. 
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I. INTRODUCTION 



Composite materials of enhanced mechanical strength, thermal conductivity or electrical 
conductivity can be produced by mixing polymer resins with nanoparticles. In particular 
carbon-based particles such as carbon nanotubes, graphene and graphite nanoplatelets are 
promising fillers because of their extraordinary materials properties l|-|3j. In this context 
percolation, a topic that had originally been brought up in the context of fluid flow through 
porous media, has attained a new field of application. 

The term percolation (more specifically connectivity percolation) refers to the transition 
at which a system changes from containing isolated clusters of particles to containing a 
system-spanning network that produces connectivity on a macroscopic scale jj, 5|. The 
density at which this transition occurs is called percolation threshold. A composite material 
made from an insulating matrix that is filled with conductive particles becomes conductive 
itself when the filler percolates 0. 

For technological applications it is often desirable to use as low a filler content as possible 
(e.g. because the filler material might be expensive, or it might have a deteriorating influence 
on other properties of a composite such as its transparency). As it has been argued that 
the percolation threshold of anisotropic objects should scale as their inverse aspect ratio Jfj], 
one commonly uses highly anisotropic fillers such as metal fibres or carbon nanotubes. For 
rod-like fillers this has been shown to be a successful approach (see e.g. jl, 8]). 

Recently plate-like fillers have moved into the focus, inparticular due to the development 
of graphite nanoplatelet reinforced polymer composites |3|, |9[ . While graphite nanoplatelets 
share some of the promising materials properties of carbon nanotubes, the production of 
graphite nanoplatelets is in general less expensive and energy intensive. However, the per- 
colation thresholds that were reached are also not as low as for carbon nanotubes. The 
authors of a recent review j^] remark that experiments do not show the expected drop of 
the percolation threshold with increasing aspect ratio. 

The topic of percolation in three dimensional (3d) suspensions of platelets has been 



addressed in theoretical and numerical investigations 



10H16I]. However, while much attention 



has been given to fully penetrable platelets, impenetrable (or otherwise interacting) platelets 
1 In contrast to the case of fluid flow through percolating cavities, however, the filler particles do not need 

to touch each other. Their distances only need to be small enough to allow for tunneling of a sufficiently 

large number of charges. 



have, to our knowledge, only been studied in two cases: Ambrosetti and co-workers did a 



simulation study on percolation in hard, oblate ellipsoids 



121 ] . and Otten and van der Schoot 



addressed the issue briefly in an article on connectedness percolation theory of polydisperse 



fillers [16|. We will refer to their results in more detail in comparison to our results. 

In this paper, we present a Monte Carlo simulation study on connectivity percolation 
in a suspension of hard platelets in 3d. We discuss, in particular, the dependence of the 
percolation threshold on the aspect ratio of the platelets, and the interference between the 
isotropic-nematic transition and the percolation transition. We interpret the simulation data 
by means of a contact-volume argument based on an effective single-particle cell model. 

II. SIMULATION DETAILS 

As model platelets we used cut spheres. These objects are obtained when a sphere of 
diameter D is intersected with two planes parallel to the equatorial plane at a distance L/2 
such that the sphere's caps are sliced off 17|. The interaction potential between two platelets 
is infinite if the platelets overlap and zero otherwise. With this interaction potential the 
system is purely entropic, its configurational properties do not depend on the temperature 
T. We use L as the unit of length and ksT as the unit of energy, where ks is Boltzmann's 
constant. 

Percolation occurs when the platelets form a system spanning cluster. A cluster is defined 
as two or more platelets which are connected. We consider two platelets as connected if 
their surfaces approach closer than a given value A. As we are using periodic boundary 
conditions, we regard a cluster as percolating if its particles are connected to their own 
periodic images. In a quantitative comparison to an experiment on electrical conductivity, 
the distance A would correspond to the tunneling distance. However, the value of A does not 
have a qualitative effect on the results we present, as long as it is smaller than the particle 
thickness L. 

Monte Carlo simulations were performed for fixed volume V and number of particles N. 
The volume of one particle is given by 

and the volume fraction is defined by 77 = NV P /V. A cuboid simulation box of dimensions 
B x B x B with periodic boundary conditions was used. Simulations for systems with a 



connecting distance A = 0.2L involved 1.5 x 10 6 trial moves per particle for equilibration 
and 3 x 10 6 moves per particle for sampling. For the other values of A equilibration and 
sampling were performed each with 0.5 x 10 6 moves per particle. 

The nematic order parameter S was used to monitor the overall alignment of the platelets 



18[. If Uj is a unit vector perpendicular to the surface of platelet i, then S is given by the 



largest eigenvalue of the tensor 



1 N 



where I is the identity matrix. In a fully aligned configuration S is unity and in a perfectly 
isotropic configuration it is zero. 

III. RESULTS AND DISCUSSION 

In an infinite system the percolation transition is marked by a discontinuity; the probabil- 
ity P(rj) to find a percolating cluster jumps from zero to one at the critical volume fraction 
T] 1 ^ . In a finite system, the transition is "smeared out", P(n) forms a sigmoidal curve, which 
becomes steeper with increasing system size (see e.g. ref. |4|). In Fig. [T]we show the perco- 
lation probability for disks with aspect ratio D/L = 5 as a function of volume fraction for 
several values of the box dimension B. The curves can be fitted by a profile of the form 

/=i|l + tanh[a(o;-6)]j , (2) 

where a is the inverse width of the transition region. In order to extract r]™ { from the 



simulation data we employ the finite-size scaling relation |4l. Il9|. 



a oc 



B 1 '" , (3) 



where v is the correlation length exponent. Fig. [2] shows the inverse width a for various 
system sizes. From the fitted line we find v = 0.96 ± 0.01, which is in agreement with the 
value for percolation in three dimensions (as well as with previous results on spheres and 



oblate ellipsoids |l2l. Il9|). 

The volume fractions at the inflection points of the tanh-profiles (eq. |2]) can be interpreted 
as finite-size percolation thresholds tj c (B). Fig. [3] shows r) c (B) plotted versus B" 1 ^. We 
fitted a line to the data using the value of v obtained above and extrapolated the fit to 
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FIG. 1. Percolation probability as a function of volume fraction r/ for different box sizes B and 
D/L = 5 (N is the number of disks inside the box). The lines are fits by a tanh-profile, the solid 
vertical line marks the infinite system percolation threshold as obtained by a finite-size scaling 
analysis. 

77™ f = 0.2419 ± 0.0001. As the difference between r] c (B) and rj 1 ^ is of order 10~ 3 , and as 
we are interested in the qualitative behaviour of the percolation threshold with the aspect 
ratio, we used the volume fraction at the inflection point as an estimate for ?7* nf for all 
other simulation data presented in the following. (Often, instead of doing a finite-size 
scaling analysis on the locations of the inflection points, one uses the intersection points of 
the percolation curves as an estimate for the percolation threshold. In Fig. [1] we marked 
r/™ f = 0.2419 by a vertical line. 77* coincides within 10~ 4 with the volume fractions of the 
intersection region.) 

In contrast to the rod- like systems which have been studied in the past, platelets ex- 
hibit an interference of the percolation transition with the isotropic-nematic (IN) transition, 
which leads to non-trivial behavior at the percolation threshold. In Fig. HI we show the 
percolation probabilty P(rj) as a function of the volume fraction rj for various aspect ratios 
D/L (connectivity distance A/L = 0.2). The volume fractions at the corresponding IN 
transitions are marked by dashed lines. For smaller aspect ratios {D/L < 8), percolation 
occurs at much smaller volume fractions than nematic ordering. For D/L = 12 the onset of 

5 
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FIG. 2. Inverse transition width a as a function of system size B for D/L = 5. The straight line 
is a fit to the data points. Its slope corresponds to a value of v = 0.96 ± 0.01. 



nematic order and the onset of percolation almost coincide. For larger D/L the situation is 
reverse: nematic ordering occurs now at much smaller volume fractions than percolation. 

Our central result is shown in Fig. [5j which depicts the percolation threshold, Tj c , as a 
function of the aspect ratio D/L for various connectivity distances A/L = 0.2...1. For 



a given A/L, r] c begins at the known value for hard spheres (D/L = 1) [20|, |2l|, then 
quickly rises to a maximum at D/L ~ 3 and upon further increasing D/L falls off weakly 
- until it reaches a minimum at intermediate aspect ratios. Upon further increasing D/L, 
t] c moderately rises again before finally reaching a plateau value. The effect of increasing 
the connectivity distance consists in shifting the curves r] c (D/L) to lower values and shifting 
the location of the minimum to higher values of D/L. The occurence of the mininum 
and the subsequent plateau-like behavior for higher D/L are consequences of the isotropic- 
nematic (IN) transition pre-empting the percolation transition. This can be seen e.g. when 
comparing Figs. |4]and|5]for D/L = 12 and A/L = 0.2. Thus when the behavior of t] c (D/L) 
beyond the minimum is analyzed one should take into account the nematic order in the 
system. In order to develop an understanding, we discuss the different regimes in detail in 
turn. 
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FIG. 3. Percolation threshold as a function of B l l u . The line is a fit to the data which yields 
rff = 0.2419 ± 0.0001. 

A. The isotropic regime 



Based on results for needle-like filler particles one would expect a decrease of r\ c with 



increasing aspect ratio 



flls 



221 ] . This kind of behavior is also predicted by a popular contact- 



volume argument which, in neglect of specific orientational correlations, assumes that r\ c is 
proportional to the ratio between the hard particle volume V and the overlap volume V^ ver 
(second virial coefficient) of two hard particles The proportionality constant a can be 
loosely associated with the number of next neighbours (contacts) within the overlap volume 
around one particle Q- 

For hard disks we have V(D, L) = V p (Eq. flQ) and V ovcv (D, L) = (%D 3 /3)(cos0(l + 
sm 2 #) + 3(cos# + [#sin#]/2)(cos#+[sin 2 #]/2) with 9 = arccos(L/£>) £]. The connectivity 
criterion set by the minimal distance A between two particles can be incorporated by taking 
V over (D + A, L + A) such that 



Vc 



2a 



V(D,L) 



(4) 



V ovei (D + A,L + A) 

For A < L < D, the contact-volume argument predicts r] c oc L/D, i.e. an inverse propor 
tionality with the aspect ratio, similar to the case of rods. 




FIG. 4. Percolation probability as a function of volume fraction for different aspect ratios (A = 
0.2L). The symbols and solid lines show the percolation probability, the dashed vertical lines 
indicate the volume fractions were the nematic order parameter S = 0.5 (approximate location of 
the isotropic- nematic transition). 

However, such a behavior is completely absent from the data. Our data for the isotropic 
phase {D/L below the minimum of i] c (D/L)) are not in the asymptotic regime in which 
Eq. (jlj) could be applied (see Fig.|5]). Using the full formula for the contact-volume argument 
(solid lines in Fig. |5]), the existence of a maximum near the sphere limit is predicted correctly, 
however, the subsequent decay with D/L is much too strong. Only for A = L one sees a 
noticeably weaker Z)/L-dependence. As a result we can state that neglecting orientational 
correlations leads to a prediction of the behavior of r] c (D/L) that is qualitatively wrong. 

These findings can be checked against the Ornstein-Zernike approach to the connected- 



ness correlation function. In a study on percolation of hard spheres [2l| the so-called "bare 
chain sum approximation" is introduced which amounts to the following overlap volume 
criterion: 

/C Kver (D + A, L + A) - Kver (A L) ' ^ ' 

It suffers from an unphysical divergence of rj c when A — > 0. For needle-like particles (sphe- 
rocylinders of length L, radius R, connectivity distance 2A with L ^> R) both the contact 
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FIG. 5. Percolation threshold as a function of disk aspect ratio for various connectivity distances. 
The solid lines represent the percolation threshold according to the contact volume ar gum ent 
(Eq. a = 1.7). The dotted straight line is the constant value of rj c predicted in Ref. [Iff] for 
the case A/L = 1. 



volume argument and the bare chain sum approximation (jSJ) give the same dependence 
r] c oc 1 / L but the volume quotient in the bare chain sum approximation is larger by a factor 
(1 + A/R)/(A/R). For disks, however, the bare chain sum approximation (D ^> L) gives 
rj c oc L/A (independent of D) instead of r/ c oc L/D from the contact volume ar gum ent. This 



has been found also in Ref. [24j and is consistent with the recent work in Ref. [16| using the 
Ornstein-Zernike approach in the isotropic fluid. Here in the limit A i C D the same 
-D-independent behavior is found: ?7 Ci oz ~ 2L/[A(5tt + 6)]. Indeed for A = L this constant 
seems to be the correct first-order approximation (see dotted line in Fig. \5§ but the other 
values for A are apparently not compatible with the domain of validity (A w L) for that 
result. 

The result of a percolation threshold that is independent of aspect ratio is in line with 
experimental observations However, as the experimental samples are in general neither 
homogeneous nor in equilibrium, it is probably not the full explanation of this effect. 
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B. The nematic regime 



For high nematic order (S > 0.9) the critical volume fraction rj c (D/L) for percolation is 
in the plateau region (see Fig. [5]). For states with high order, one can invoke a single-particle 
cell model to derive the following expression which connects volume fraction 77, aspect ratio 
D/L and order parameter S 

7 D 



l + P-VT^S . (6) 
77 L 

In the cell model, one assumes that each disk occupies an effective volume V e s = (it / 'A)D 2 L e g 
and the effective height of the disk is calculated via the one-particle orientation distribution 



function (see App. IX]) . L cS = L + 0Dy/l — S. We determined the constant /3 as 0.49 by 
a fit to our data (see Fig. E]), being somewhat smaller than derived from simple guesses 
for the orientation distribution (/3 ~ 0.7). The constant 7 = 0.88 is close to 1, which is 
consistent with the interpretation that the effective volume for all disks almost fills space. A 



cell model description is known to be fairly accurate for the columnar phase of disks [25|, [26 ] 
where one can go further and also derive the full partition function. For the percolation 
problem, we are in the nematic phase but still far away from the columnar phase. Hence it 
is somewhat surprising that an effective one-particle argument quantifies the value of the 
order parameter for given volume fraction and aspect ratio also in the nematic phase. 

In order to explain the plateau in Fig. \5[ we make the assumption that the onset of 
percolation is dictated by A oc L c r — L. It is a similar argument as the contact volume 
argument in the isotropic phase, namely that there is an invariant for percolation which is a 
quotient between two volumina: the effective excess volume of one disk (7r/4).D 2 (L e g- — L) and 
the excess overlap volume of two aligned disks ~ (tt/A)D 2 A. Using Eq. ([6]), we immediately 
find for the critical volume fraction at percolation: 

1 = 1 + P , T • (7) 

With this assumption, t] c is indeed independent of D in the nematic phase, as seen in the 
data. In Fig. [6] we plot the plateau values of l/r] c vs. A/L and also observe the predicted lin- 
ear dependence. A linear fit, however, gives 7 ~ 0.61 (in contrast to 7 ~ 1 from the equation 
of state (E]). The value of 7 is the critical volume fraction for A = and thus corresponds 
to the volume fraction for the jammed state. This simple extrapolation, however, neglects 
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(1-S) 1/2 D/L A/L 

FIG. 6. Left panel: Plot of inverse volume fraction I/77 vs. the scaling variable y/1 — SD/L for 
our data deep in the nematic phase with order parameter S > 0.9. The linear fit assumes the 
cell theory equation Q. Right panel: Plot of inverse critical volume fraction l/r] c at the plateau 
(D/L S> 1) vs. connectivity parameter A/L. The linear fit is according to the volume argument 
in Eq. ©. 

the fact that the system undergoes a transition to the columnar state before reaching such 
high values for the volume fraction. 

When the system is in the isotropic phase, the percolating cluster is connected to its 
periodic images in all spatial directions with equal probability. As the system undergoes 
the IN transition an asymmetry develops. In Fig. [7] we show the ratio of clusters that 
percolate "along the director" (in the z-direction) to those that percolate "perpendicular" 
to the director, i.e. that are connected to their periodic images in the x- and/or ^/-direction. 
Once the system enters nematic phase, percolation along the director becomes less likely. 



IV. SUMMARY AND CONCLUSIONS 

We have studied connectivity percolation for systems of monodisperse hard platelets 
modeled by cut spheres of diameter D and height L. The connectivity criterion between 
two particles was modeled by a minimum surface-to-surface distance A chosen to be of 
the order of the particle height, L ~ A. For a fixed A, we have found that the critical 
volume fraction r\ c for percolation varies surprisingly little with the aspect ratio D/L (see 
Fig. [S]). The most notable feature of rj c (D) is a shallow local minimum followed by a plateau 

11 



I 

J 

|D/L=15 
1.5 + 



+ 



1 J — •••\V . . Jgfe^ta*** * ..[.......*..• 

- 1 A A A A 1 

3TA A A A | A- 

I I 



• D/L = 5 

a D/L = 15 



).2 



0.3 



0.4 



0.5 



0.6 



0.7 



FIG. 7. Ratio of clusters that percolate along the director to clusters that percolate perpendicular 
to it. The vertical lines mark the IN transition (S(rj) = 0.5). 



for D/L ^> 1. We have demonstrated that the appearance of the minimum is directly 
connected to the onset of nematic ordering in the system. The Z)-independent, plateau- 
like behaviour of r\ c and the variation of the value of r\ c on the plateau with A could be 
rationalized in terms of a contact-volume argument (see Sec. IIII B[) . This argument relies 
on a derivation of the dependence of the nematic order parameter S on aspect ratio D/L 
and volume fraction 77, based on an effective single-particle cell model. We have tested this 
relation by our simulation results and found good agreement. In contrast to the nematic 
phase, a simple interpretation of rj c (D/L) in the isotropic phase (small aspect ratios) cannot 
be given, as angular correlations appear to play a major role. A theoretical approach based 
on the solution of integral equations for the connectivity correlations appears to be promi sing 



in this respect 16], further work based on fundamental measure density functionals 271. |28| 



and related reference system closures 



29 



30] is desirable in our opinion. 



It is a common strategy to use highly anisotropic, conductive particles as fillers in com- 
posite materials to achieve conductivity at low filler concentrations. This strategy is based 
on a "rule of thumb" which says that t] c decreases with the aspect ratio of anisotropic par- 
ticles. A major conclusion from our work is that this rule does not hold for platelets such 



as graphite nano-platelets. This seems to be in contradiction to the results in Ref. 



12| 



12 



where percolation of hard ellipsoids is studied and the corresponding rj c for ob 



(i.e. disks) shows a strong decrease with aspect ratio (see Fig. 3 in Ref. 12(). However 



ate ellipsoids 



these data are obtained as curves rj c (D/L) for constant A/D (and not A/L as in our case), 
i.e. along these curves the connectivity distance increases with aspect ratio. We have re- 
analyzed the data in terms of a family of curves for constant A/L and found that apparently 
all data have been obtained in the isotropic phase and there is little variation of rj c (D/L), 
similar to our findings in the isotropic phase. Hence the authors of Ref. [3] could not see 
the interference of the IN transition with the percolation transition. 
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Appendix A: Cell model for strong nematic order 

We consider disks deep in the nematic state, i.e. S ^ 0.9, see Fig. [H] for geometric 
definitions for one disk. The deviation of the symmetry axis of one disk (disk orientation) 
with respect to the director is given by 9. The order parameter S is defined by 

s /W0-1\ (Al) 



where (...) denotes the thermal average. We call p(9) the probability distribution function 
of disk orientations in the bulk state. Using this, the order parameter becomes 

1 f , m 3cos 2 0-l 1 f 1 „ „, .„,3cos 2 0-l j . v 

S = 4W dnp{6) 2 = 2 / x (C ° S )P{ ) 2 • (A2) 

Through fluctuations around the director, each disk will occupy approximately an effective 
volume (D ^> L): 

V cS = ^D 2 L cS = ^D 2 (L + Dsm9) . (A3) 
13 



The assumption of our cell model is that the whole system can be considered to be an 
ensemble of closely packed effective volumes which are occupied by single disks which move 
independently of each other. In order to determine that effective volume we consider a 
parametrization of the probability distribution function in the form p(9) = N a p Te ^(9/a) 
where a is a width parameter in the angular variable and N a a normalization constant. 
Furthermore we consider narrow distributions (large a) which permit us to approximate 
angular averages as 

i r 1 

2 

where the dots stand for any angular variable. Let us define moments of the reduced dis- 
tribution function by M n = J °° x n p ie ^(x)dx. Using these definitions, we have the following 
relations: 



/l poo 
d(cos9)...p(9) « J 9d9...p(9), (A4) 



D 



poo 

1 = / 9d9p{9) = N a a 2 M x (normalization) , (A5) 
Jo 

T J POO 

- / 9 2 d9p(9) = N a a 3 M 2 (effective height) , (A6) 
Jo 



OO q 

3ja„ra\ - 6 AT „A 



S = - 9 6 d9p{9) = -N a a 4 M 3 (order parameter) . (A7) 
2 ./ 2 



We can use the equations for normalization and order parameter to eliminate a and and 
find 

L cS = L + D PVT^S (A8) 

Mi V 3 Mi 

where (3 depends only on the type of reduced distribution p re d- As examples, for an expo- 
nential distribution p rc d{x) = exp(— x) (as found in the cell theory for the columnar phase 



25 



261 ]) we have (3 = 2/3 and for a Gaussian distribution p rc d{x) = exp(— x 2 ) we find 
(3= v /^76^0.72. 

As stated before, we proceed with the assumption that the sum of the effective volumes 
V e s for particles (almost) fills the system volume V, i.e. NV c g/V — 7 ~ 1. Using the 
volume fraction rj = N(n / 4)D 2 L/V (D 3> L) we rewrite this as r)V e s/[(ir/4)D 2 L] = 7 and 



upon using Eqs. (IA3j) and ( 1A8j) we find the final result 

7 = 7/(l + /3v / I^s£) • (A9) 



14 



DsinG 



L 



D 



FIG. 8. Geometrical definitions. 



[1] M. Moniruzzaman and K. I. Winev. iMacromoleculesI 39. 5194 (2006), ISSN 0024-9297 

[2] S. Stankovich, D. Dikin, G. Dommett, K. Kohlhaas, E. Zimney, E. Stach, R. Piner, S. Nguyen, 

and R. Ruoff, Nature 442, 282 (2006), ISSN 0028-0836 
[3] B. Li and W. Zhong, J. Mat. Sci. 46, 5595 (2011), 

http : //www. springerlink. com/ content/g4p28p38t285830v/ 



[4] D. Stauffer and A. Aharony, Introduction To Percolation Theory, 2nd ed. (Taylor & Francis, 

London, 1994) |http : //amazon . com/o/ASIN/B000PUB6Wq/| 
[5] R. Meester and R. Roy, Continuum Percolations (Cambridge University Press, Cambridge, 

1996) 

[6] I. Balberg, C. Anderson, S. Alexander, and N. Wagner, Phys. Rev. B 30, 3933 (1984), ISSN 
1550-235X, |http : //prb . aps . org/abstract/PRB/v30/i7/p3933_l 



[7] H. Deng, R. Zhang, E. Bilotti, J. Loos, and T. Peijs, J. Appl. Polym. Sci. 113, 742 (2009) 

[8] K. Leung and D. Chandler, J. Stat. Phy. 63, 837 (1991) 

[9] J. Potts, D. Dreyer, C. Bielawski, and R. Ruoff, Polymer 5, 25 (2011) 

[10] E. Charlaix, J. Phys. A: Math. Gen. 19, L533 (1986) 



[11] Y. B. Yi and E. Tawerghi, Phys. Rev. E 79, 041134 (Apr 2009) 



|http : //pre . aps . org/abstract/PRE/v79/i4/e041134 



[12] G. Ambrosetti, N. Johner, C. Grimaldi, A. Danani, and P. Ryser, Phys. Rev. E 78, 61126 



(2008), ISSN 1550-2376, http : //pre . aps . org/abstract/PRE/v78/i6/e061126 



[13] J. Quintanilla and R. Ziff, Phys. Rev. E 76, 51115 (2007), ISSN 1550-2376, 
|http : //pre . aps . org/abstract/PRE/v76/i5/e051115 



15 



[14] J. Quintanilla, Phys. Rev. E 63, 61108 (2001), ISSN 1550-2376, 

|http : //pre . aps . org/abstract/PRE/v63/i6/e061108| 

[15] J. Quintanilla, S. Torquato, and R. Ziff, J. Phy. A: Math. Gen. 33, L399 (2000), 

|http : //iopscience . iop . org/0305-4470/33/42/104 



[16] R. Otten and P. van der Schoot, J. Chem. Phys. 134, 094902 (2011) 



http : // j cp . aip . org/resource/1/ j cpsa6/vl34/i9/p094902_sl 



[17] J. Veerman and D. Frenkel, Phys. Rev. A 45, 5632 (1992), ISSN 1094-1622, 
|http : //pra. aps . org/abstract/PRA/v45/i8/p5632_l 



[18] P. G. de Gennes and J. Prost, The physics of liquid crystals (Oxford University Press, Oxford, 
1993) 

[19] M. Rintoul and S. Torquato, |J. Phy. A: Math. Gen.|30, L585 (1997) 



[20] A. Bug, S. Safran, G. Grest, and I. Webman, Phys. Rev. Lett. 55, 1896 (1985), ISSN 1079- 



7114, http : //prl . aps . org/abstract/PRL/v55/il8/pl896_l 



[21] T. DeSimone, S. Demoulini, and R. Stratt, J. Chem. Phys. 85, 391 (1986), 



http : // j cp . aip . org/resource/ 1/ j cpsa6/v85/ i l/p391_sl 



[22] T. Schilling, S. Jungblut, and M. Miller, Phys. Rev. Lett. 98, 108303 (2007), ISSN 1079-7114, 
|http IT/ prl . aps . org/ abstract/PRL/v98/ il0/el08303 



[23] D. Frenkel, Liq. Cryst. 5, 929 (1989) 

[24] G. Ambrosetti, C. Grimaldi, I. Balberg, T. Maeder, A. Danani, and P. Ryser, Phys. Rev. B 
81, 155434 (2010), |http : //prb . aps . org/abstract/PRB/v81/il5/el55434 



[25] H. Wensink, Phys. Rev. Lett. 93, 157801 (2004) 

[26] H. Wensink and H. Lekkerkerker, Mol. Phys. 107, 2111 (2009) 

[27] A. Esztermann, H. Reich, and M. Schmidt, Phys. Rev. E 73, 011409 (2006), 
|http: //pre . aps . org/abstract/PRE/v73/il/e0ii409 



[28] J. Phillips and M. Schmidt, Phys. RevTE] 81, 041401 (Apr 2010) 



http : //pre . aps . org/abstract/PRE/v81/i4/e041401 



[29] M. Oettel, J. Phys.: Condens. Matter 17, 429 (2005) 

[30] A. Ayadim, M. Oettel, and S. Amokrane, J. Phys.: Condens. Matter 21, 115103 (2009) 



16 



